# Required packages
require(tidyverse)
require(ggthemes)
require(flextable)

# Load data and models
load("01_Data/spatial_dimension_models.Rdata")
load("01_Data/spatial_dimension_data.Rdata")
load("01_Data/spatial_dimension_weight_matrix.Rdata")

#### Figure A1: Weight matrix visualized ####
jpeg("02_Figures/appendix_A1.jpg", width = 1000, height = 1200)
par(bg = "white")
plot(invd.weights,
     location_coordinates,
     lwd = .05,
     col = "black",
     cex = .1)
dev.off()


# Source helper functions
# These contain the calculation of the different standard errors 
# and scripts for the creation of the final tables

source("00_Scripts/00_helper_function.R")

#### Table 1 and Table A2: Main regression results ----
models <- list(lm1, lm2, lm3,
               sm1, sm2, sm3,
               lmer1 , lmer2, lmer3)

models_df <- bind_rows(lapply(models,get_model_df))

get_table(models_df, file_name = "03_Tables/table1.docx", full_table =F)
get_table(models_df, file_name="03_Tables/table1_full.docx", full_table =T)

#### Figure 2: Coefficient plot ----
pdf("02_Figures/Figure02.pdf", height=5, width=10)
models_df %>%
  filter(term == "Distance (hours)" |
           term == "Risk (high)" |
           term == "Risk (medium)") %>%
  filter(id != "3: Linear (OLS)" &
           id != "6: SAR (MLE)"  &
           id != "9: LME (MLE)") %>%
  ggplot(aes(
    x = estimate,
    y = term,
    color = id,
    shape = id,
    xmax = conf.high,
    xmin = conf.low
  )) +
  geom_linerange(aes(xmax= conf.9.high,
                     xmin = conf.9.low),
                 position = position_dodge2(width = .4,reverse = T),
                 size =.9) +
  geom_pointrange(position =  position_dodge2(width = .4, reverse=T),
                  size = .5) +
  xlab("Coefficient") +
  ylab("") +
  ggthemes::theme_base() +
  scale_color_brewer(palette = "Dark2") +
  geom_vline(xintercept = 0, linetype=2) +
  theme(legend.title = element_blank(),
        plot.background =  element_blank())
dev.off()

#### Figure A3.1 - A3.3: Random effect plot ----
pdf("02_Figures/appendix_A3_1.pdf", width=18, height=7)
broom.mixed::tidy(lmer1, "ran_vals", conf.int = TRUE) %>% 
  filter(group == "qtr:district") %>% 
  separate(level, c("qtr","district"),":") %>% 
  ggplot(aes(x=qtr, color=district, y=estimate, ymax=conf.high, ymin=conf.low))+
  geom_hline(yintercept =  0, linetype = 2)+
  geom_pointrange(alpha=1, position=position_dodge(width=.5))+
  xlab("")+
  ylab("Random effect (Intercept)")+
  theme_base()+
  scale_color_brewer(palette = "Dark2") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position="bottom",
        legend.title = element_blank(),
        plot.background = element_blank())
dev.off()



pdf("02_Figures/appendix_A3_2.pdf", width=18, height=7)
broom.mixed::tidy(lmer2, "ran_vals", conf.int = TRUE) %>% 
  filter(group == "qtr:district") %>% 
  separate(level, c("qtr","district"),":") %>% 
  ggplot(aes(x=qtr, color=district, y=estimate, ymax=conf.high, ymin=conf.low))+
  geom_hline(yintercept =  0, linetype = 2)+
  geom_pointrange(alpha=1, position=position_dodge(width=.5))+
  xlab("")+
  ylab("Random effect (Intercept)")+
  theme_base()+
  scale_color_brewer(palette = "Dark2") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position="bottom",
        legend.title = element_blank(),
        plot.background = element_blank())
dev.off()

pdf("02_Figures/apendix_A3_3.pdf", width=18, height=7)
broom.mixed::tidy(lmer3, "ran_vals", conf.int = TRUE) %>% 
  filter(group == "qtr:district") %>% 
  separate(level, c("qtr","district"),":") %>% 
  ggplot(aes(x=qtr, color=district, y=estimate, ymax=conf.high, ymin=conf.low))+
  geom_hline(yintercept =  0, linetype = 2)+
  geom_pointrange(alpha=1, position=position_dodge(width=.5))+
  xlab("")+
  ylab("Random effect (Intercept)")+
  theme_base()+
  scale_color_brewer(palette = "Dark2") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position="bottom",
        legend.title = element_blank(),
        plot.background = element_blank())
dev.off()


rm(list = ls())
